DF检验

考虑如下的基础模型 \(x_t = \phi_0+\phi_1x_{t-1}+\xi_t, \xi_t\)为随机部分，且\(\xi_t\sim N(0,\sigma^2)\)

提出假设：原假设为序列非平稳\(H_0:\vert\phi_1\vert\ge1\),备择假设是序列平稳\(H_1:\vert\phi_1\vert\lt1\)

检验统计量： \[DF=\frac{\hat\phi_1-1}{SE(\hat\phi_1)}\] \(DF\)检验模型有三种： - 无漂移项自回归模型：\(x_t=\phi_1x_{t-1}+\xi_t\) - 有漂移项自回归模型：\(x_t = \phi_0+\phi_1x_{t-1}+\xi_t\) - 带趋势回归模型：\(x_t = \alpha + \beta t+\phi_1x_{t-1}+\xi_t\)

ADFTest

考虑如下的基础模型： \[X_t = c_t+\beta X_{t-1}+\sum_{j=1}^{p-1}\phi_j\Delta X_{t-j}+e_t\] 当\(\beta=1\)时，就是\(\Delta X_t\)的AR(p-1)模型，当\(\beta<1\)时，是\(X_t\)的AR(p)模型。 \(c_t\)是非随机的趋势部分，可以取0，或常数，或a+bt这样的非随机线性趋势。

检验假设: \(H_0:\beta=1\) vs. \(H_1:\beta<1\)

如果拒绝H0，就说明没有单位根，使用统计量\[ADF=\frac{\hat\beta-1}{SE(\hat\beta)}\] 当ADF统计量足够小的时候拒绝H0.

• fUnitRoots包的adfTest()函数可以执行单位根ADF检验。
• tseries包的adf.test()函数也可以执行单位根ADF检验。

注意，单位根DF检验和ADF检验都是在拒绝（显著）时否认有单位根， 不显著时承认有单位根。

PP检验：

ADF和DF检验要求随机扰动项独立同分布或者独立同正态分布，但是很多时间序列不满足这一点，故此提出\(PP\)检验，它限定随机序列可以具有暂时相关性和不同分布，它的原假设时间序列含有一个单位根，对于模型\(x_t = \phi x_{t-1}+\xi_t\). 提出原假设：\(H_0:\phi=1\)

ADF检验主要适用于方差齐性的场合，它对于异方差序列的平稳性检验结果不佳。

PP检验的适用于异方差场合的单位根检验，且PP检验的临界值和ADF检验相同,服从相应的ADF检验统计量的极限分布。

KPSS检验

KPSS检验：ADF检验和PP检验的原假设都是序列中有单位根，也就是认为序列是不平稳的，备择假设是序列是平稳的。在实证研究中，人们愿意多做差分（有单位根就需要做差分）而不愿意忽略单位根的存在，因为忽略单位根比过度差分带来的后果更为严重。其次以单位根作为原假设相对更加容易涉及检验统计量。故此以序列的平稳性为原假设的方法产生：KPSS方法。

单位根检验实例：

da <- read_table2("~/Downloads/ftsdata/q-gnp4710.txt", col_types=cols(
  .default = col_double()))
gnp <- ts(log(da[["VALUE"]]), 
          start=c(1947, 1), frequency=4)
dgnp <- diff(gnp)
rm(da)
plot(dgnp, xlab="year", ylab="log(GNP)")

forecast::Acf(dgnp, main="")

forecast::Pacf(dgnp, main="")

#pacf截尾，acf拖尾，考虑用ar模型建模
#用aic定阶
ar(dgnp,method='mle')

## 
## Call:
## ar(x = dgnp, method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.4318   0.1985  -0.1180   0.0189  -0.1607   0.0900   0.0615  -0.0814  
##       9  
##  0.1940  
## 
## Order selected 9  sigma^2 estimated as  8.918e-05

#AIC取9阶
#ADF的基础模型需要一个AR阶数,取p=9。 
#用fUnitRoots::adfTest()对GNP的对数值进行ADF单位根检验：
library(fUnitRoots)
#type: c: with constant(intercept)
#     nc: no constant
#     ct:constant and trend
#这里gnp的均值非0，且存在trend，故type='ct'
fUnitRoots::adfTest(gnp, lags=9, type="ct") 

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 9
##   STATISTIC:
##     Dickey-Fuller: -0.0094
##   P VALUE:
##     0.99 
## 
## Description:
##  Tue Mar  9 00:19:35 2021 by user:

fUnitRoots::unitrootTest(gnp, lags=9, type="ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 9
##   STATISTIC:
##     DF: -0.0094
##   P VALUE:
##     t: 0.996 
##     n: 0.996 
## 
## Description:
##  Tue Mar  9 00:19:35 2021 by user:

#结果值较大，说明不能拒绝零假设， 即对数GNP序列有单位根。序列非平稳
#另外两种单位根检验方法
#tseries::pp.test(gnp)
#tseries::kpss.test(gnp)
#尝试人为地拟合非随机线性增长趋势，检验残差是否有单位根
tmp.t <- c(time(gnp))
tmp.y <- residuals( lm(c(gnp) ~ tmp.t) )
fUnitRoots::adfTest(tmp.y, type="nc")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.0763
##   P VALUE:
##     0.592 
## 
## Description:
##  Tue Mar  9 00:19:35 2021 by user:

#结果说明用回归去掉非随机的线性增长趋势后仍有单位根存在。

SARIMA

da <- read_excel('~/Desktop/应用时间序列实验课/实验课/数据/中国GDP数据.xlsx',sheet=2)
da <- ts(da$`季度GDP(亿元)`,frequency = 4,start=c(2008,1))
autoplot(da)

ggseasonplot(da)

forecast::ggsubseriesplot(da)

ggtsdisplay(da)

#趋势，周期
da %>% diff(1) %>% ggtsdisplay()

da %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

da %>% diff(lag=12) %>% diff() %>% ggtsdisplay()

自动拟合arima的阶数

arima.model <- auto.arima(da,approximation = F);arima.model

## Series: da 
## ARIMA(0,1,2)(1,1,0)[4] 
## 
## Coefficients:
##           ma1     ma2     sar1
##       -0.8613  0.2006  -0.4118
## s.e.   0.1209  0.1261   0.1140
## 
## sigma^2 estimated as 33332300:  log likelihood=-694.77
## AIC=1397.54   AICc=1398.17   BIC=1406.48

#SARIMA(p,d,q)(P,D,Q)[s]
#S:季节周期的长度
#D:季节差分的阶数

\[(1+0.4118B^4)(1-B)(1-B^4)y_t = (1-0.8613B+0.2006B^2)\epsilon_t\]

resm <- Arima(
  da, order=c(0,1,2), seasonal=list(order=c(1,1,0), period=4)
); 
summary(resm)

## Series: da 
## ARIMA(0,1,2)(1,1,0)[4] 
## 
## Coefficients:
##           ma1     ma2     sar1
##       -0.8613  0.2006  -0.4118
## s.e.   0.1209  0.1261   0.1140
## 
## sigma^2 estimated as 33332300:  log likelihood=-694.77
## AIC=1397.54   AICc=1398.17   BIC=1406.48
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 580.3015 5452.413 2986.756 1.229971 3.998121 0.3203815 -0.01847175

#残差检验
autoplot(resm$residuals) #非白噪声，波动聚集

Box.test(resm$residuals,lag=12,type = 'Ljung')

## 
##  Box-Ljung test
## 
## data:  resm$residuals
## X-squared = 5.7441, df = 12, p-value = 0.9284

Box.test(resm$residuals,lag=24,type = 'Ljung')

## 
##  Box-Ljung test
## 
## data:  resm$residuals
## X-squared = 11.725, df = 24, p-value = 0.9828

Box.test(resm$residuals^2,lag=6,type = 'Ljung') #非白噪声，进一步进行异方差建模

## 
##  Box-Ljung test
## 
## data:  resm$residuals^2
## X-squared = 13.657, df = 6, p-value = 0.03371

library(fGarch)
mod1 <- garchFit( ~ garch(1,1), data=resm$residuals, trace=FALSE)
#其中1表示均值方程是一个常数，
#输出结果中mu表示均值方程的均值，omega表示alpha0,alpha1为alpha1
resi <- residuals(mod1, standardize=TRUE)
checkresiduals(resi)

#标准化残差的时序图
plot(ts(resi, start=start(da), frequency=frequency(da)), 
     xlab="年", ylab="标准化残差")

#标准化残差的白噪声检验
Box.test(resi,lag=12)

## 
##  Box-Pierce test
## 
## data:  resi
## X-squared = 12.588, df = 12, p-value = 0.3997

Box.test(resi^2,lag=12)

## 
##  Box-Pierce test
## 
## data:  resi^2
## X-squared = 7.3936, df = 12, p-value = 0.8305

ggAcf(resi)

ggPacf(resi)

• 超前多步预报

tmp.y <- window(da, start=start(da), end=c(2024,4))
resm2 <- arima(tmp.y, order=c(0,1,2), seasonal=list(order=c(1,1,0), period=4))
#样本外预测
resm2 %>% forecast(h=12) %>% autoplot()

mode2 <- garchFit( ~ garch(1,1), data=resm2$residuals, trace=FALSE)
pred1 <- predict(resm2, n.ahead=6)
mod <- predict(mode2,n.ahead=6)

• 样本内静态预测（预测样本内最后6个观测）

horiz=6
dbnp.prev=predict(resm2,n.ahead=horiz)$pred
dbnp.se=predict(resm2,n.ahead=horiz)$se
bornesup=dbnp.prev+1.96*dbnp.se
borneinf=dbnp.prev-1.96*dbnp.se
library(ggplot2)
ts.plot(da)
lines(dbnp.prev,col="red")
lines(bornesup,col="orange")
lines(borneinf,col="orange")

• 样本内动态预测（预测一个前进一个建模）

dbnp.prev.dyn=ts(start=length(da)-horiz+1,end=length(da))
dbnp.se.dyn=ts(start=length(da)-horiz+1,end=length(da))
for (i in (length(da)-horiz+1):length(da)){
  
  dbnp.mod.temp=arima(window(c(da),end=i),order=c(0,1,2),seasonal=list(order=c(1,1,0),period=4),include.mean=F)
  dbnp.prev.dyn[i-length(da)+horiz]=predict(dbnp.mod.temp,n.ahead=1)$pred
  dbnp.se.dyn[i-length(da)+horiz]=predict(dbnp.mod.temp, n.ahead=1)$se
}
par(mfrow=c(1,1))
ts.plot(window(c(da),start=1,end=74),col="orange",main='Dynamic Prediction')
lines(dbnp.prev.dyn,col="red")
lines(dbnp.prev.dyn+1.96*dbnp.se.dyn,col="blue")
lines(dbnp.prev.dyn-1.96*dbnp.se.dyn,col="blue")

rmse=sqrt(mean((c(da)[(length(da)-horiz+1):length(da)]-dbnp.prev.dyn)**2))
rmse

## [1] 18377.35

ARIMA拟合

arima.model2 <- auto.arima(da,approximation = F,seasonal = F,include.mean=T);arima.model2

## Series: da 
## ARIMA(4,1,1) with drift 
## 
## Coefficients:
##           ar1     ar2     ar3     ar4      ma1      drift
##       -0.0568  0.0093  0.0422  0.8713  -0.7753  1962.7763
## s.e.   0.0629  0.0571  0.0588  0.0538   0.0892   926.1329
## 
## sigma^2 estimated as 38114765:  log likelihood=-740.57
## AIC=1495.13   AICc=1496.86   BIC=1511.17

#注：The parameter muis called the “drift” in the R output when d=1
checkresiduals(arima.model2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,1) with drift
## Q* = 15.543, df = 3, p-value = 0.001407
## 
## Model df: 6.   Total lags used: 9

Box.test(arima.model2$residuals,lag=12) #白噪声

## 
##  Box-Pierce test
## 
## data:  arima.model2$residuals
## X-squared = 16.203, df = 12, p-value = 0.1821

Box.test(arima.model2$residuals^2,lag=12) #白噪声

## 
##  Box-Pierce test
## 
## data:  arima.model2$residuals^2
## X-squared = 17.662, df = 12, p-value = 0.1263

\[(1+0.0568B-0.0093B^2-0.0422B^3-0.8714B^4)(1-B)^1(y_t-1962t)=(1-0.7753B)e_t\]

非季节的ARIMA模型可以写为： \[(1-\phi_1B-...-\phi_pB^p)(1-B)^dy_t=c+(1+\theta_1B+...+\theta_qB)e_t\] 上式表达等价于: \[(1-\phi_1B-...-\phi_pB^p)(1-B)^d(y_t-\mu t^d/d!)=(1+\theta_1B+...+\theta_qB)e_t\] 其中B是滞后算子，\(c=\mu (1-\phi_1-...-\phi_p)\),及\(\mu\)是\((1-B)^dy_t\)的均值。 当d=0,且include.drift=T的时候，模型形式为： \[(1-\phi_1B-...-\phi_pB^p)(y_t-a-bt)=(1+\theta_1B+...+\theta_qB)e_t\] a:intercept;b:drift

• SARIMA 包含了季节项的滞后算子，如对于\(ARIMA(1,1,1)(1,1,1)_4\)模型， 模型可以写作： \[(1-\phi_1B)(1-\Phi_1B^4)(1-B)(1-B^4)y_t=(1+\theta_1B)(1+\Theta_1B^4)\epsilon_t\]

相关连接：https://robjhyndman.com/hyndsight/arimaconstants/